In [1]:
%reload_ext autotime
import rasterio
import requests
from rasterio.features import shapes, sieve, rasterize
from rasterio.plot import show
import geopandas as gpd
import os
import numpy as np
from tqdm.auto import tqdm
from shapely.geometry import box, shape
from rasterio.mask import mask
from matplotlib import pyplot as plt
import xdem
import pandas as pd
tqdm.pandas()
pd.set_option('display.max_columns', None)

✔️ 2.09 s (2024-12-02T11:03:52/2024-12-02T11:03:54)
In [2]:
manifest = pd.json_normalize(requests.get("https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2023/dem_1m/2193/collection.json").json()["links"])
manifest = manifest[manifest["rel"] == "item"]
manifest["tilename"] = manifest.href.str.replace(".json", "").str.strip("./")
manifest

✔️ 211 ms (2024-12-02T11:03:54/2024-12-02T11:03:54)
Out[2]:
rel href type file:checksum tilename
2 item ./BD43_10000_0204.json application/json 12208b52bfdff84d924430e9503bc70bcc7555c1861e5b... BD43_10000_0204
3 item ./BD43_10000_0205.json application/json 12204ec4b1257e965c8da142c4696b1f8cbbc6c78d0cf0... BD43_10000_0205
4 item ./BD43_10000_0304.json application/json 1220acfd4fd86a3866f4facf92bb75cd5224451e774ff4... BD43_10000_0304
5 item ./BD43_10000_0305.json application/json 122051507adf6a042f09c138f1842911fec7bd43c270d8... BD43_10000_0305
6 item ./BD43_10000_0404.json application/json 1220c63b7011af35652df009f8a6b69f4e3f8ea2db33d6... BD43_10000_0404
... ... ... ... ... ...
302 item ./BH43_10000_0102.json application/json 12200f2d3649a677e745f258162ca6ebc13904d6ecde95... BH43_10000_0102
303 item ./BH43_10000_0201.json application/json 1220d28f9ff3b2f6ad76139b4ccfa5bd86bc996b417e36... BH43_10000_0201
304 item ./BH43_10000_0202.json application/json 1220d78d97c89b75b5d35667f77acb8e5a2078d52b977b... BH43_10000_0202
305 item ./BH43_10000_0301.json application/json 12200b058bab70e426eccd9ae990a5888e1754b27b97c2... BH43_10000_0301
306 item ./BH43_10000_0302.json application/json 12203c93b413b89941df279575d9d352046cc8b6977deb... BH43_10000_0302

305 rows × 5 columns

In [3]:
tilename = "BD44_10000_0503"
old = rasterio.open(f"https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2018-2020/dem_1m/2193/{tilename}.tiff")
new = rasterio.open(f"https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2023/dem_1m/2193/{tilename}.tiff")

✔️ 929 ms (2024-12-02T11:03:55/2024-12-02T11:03:55)
In [4]:
new.crs, old.crs

✔️ 2.56 ms (2024-12-02T11:03:56/2024-12-02T11:03:56)
Out[4]:
(CRS.from_wkt('PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","2193"]]'),
 CRS.from_wkt('PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","2193"]]'))
In [5]:
show(old, title="2018-2020 DEM")

✔️ 4.8 s (2024-12-02T11:03:56/2024-12-02T11:04:01)
No description has been provided for this image
Out[5]:
<Axes: title={'center': '2018-2020 DEM'}>
In [6]:
show(new, title="2023 DEM")

✔️ 4.88 s (2024-12-02T11:04:01/2024-12-02T11:04:06)
No description has been provided for this image
Out[6]:
<Axes: title={'center': '2023 DEM'}>
In [7]:
diff = new.read(1) - old.read(1)
print(diff.shape)
plt.imshow(diff, cmap="coolwarm_r", vmin=-1, vmax=1)
plt.colorbar()

✔️ 1.4 s (2024-12-02T11:04:06/2024-12-02T11:04:07)
(7200, 4800)
Out[7]:
<matplotlib.colorbar.Colorbar at 0x7f75e42eb430>
No description has been provided for this image
In [8]:
_ = plt.hist(diff.flatten(), bins=100, range=(-2, 2))
plt.title("Histogram of elevation differences")

✔️ 448 ms (2024-12-02T11:04:07/2024-12-02T11:04:08)
Out[8]:
Text(0.5, 1.0, 'Histogram of elevation differences')
No description has been provided for this image
In [9]:
result = diff.round().clip(min=-1, max=1).astype(np.int16)
result = sieve(result, 4000)
result = np.where(result != 0, result, np.nan)
plt.imshow(result, cmap="coolwarm_r")
plt.colorbar()

✔️ 3.23 s (2024-12-02T11:04:08/2024-12-02T11:04:11)
Out[9]:
<matplotlib.colorbar.Colorbar at 0x7f75e413f820>
No description has been provided for this image
In [10]:
areas = gpd.GeoDataFrame(geometry=[shape(s) for s, v in shapes(result, transform=new.transform) if not np.isnan(v)])
areas["area"] = areas.area
areas.sort_values("area", ascending=False, inplace=True)
areas.plot("area", legend=True)

✔️ 1.35 s (2024-12-02T11:04:11/2024-12-02T11:04:12)
Out[10]:
<Axes: >
No description has been provided for this image
In [11]:
largest_area = areas.iloc[0]#.head(1)
largest_area

✔️ 6.1 ms (2024-12-02T11:04:12/2024-12-02T11:04:12)
Out[11]:
geometry    POLYGON ((2054419 5803910, 2054419 5803909, 20...
area                                                 530308.0
Name: 176, dtype: object
In [12]:
masked_old, transform = mask(old, [largest_area.geometry], nodata=np.nan)
nonan = np.argwhere(~np.isnan(masked_old[0]))
top_left = nonan.min(axis=0)
bottom_right = nonan.max(axis=0)
print(top_left, bottom_right)
masked_old = masked_old[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_old)

✔️ 1.25 s (2024-12-02T11:04:13/2024-12-02T11:04:14)
[5290    0] [7044 3057]
No description has been provided for this image
Out[12]:
<Axes: >
In [13]:
masked_new, transform = mask(new, [largest_area.geometry], nodata=np.nan)
masked_new = masked_new[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_new)

✔️ 1.02 s (2024-12-02T11:04:14/2024-12-02T11:04:15)
No description has been provided for this image
Out[13]:
<Axes: >
In [14]:
features = {
    "old_min": np.nanmin(masked_old),
    "old_max": np.nanmax(masked_old),
    "old_mean": np.nanmean(masked_old),
    "old_median": np.nanmedian(masked_old),
    "old_std": np.nanstd(masked_old),
    "new_min": np.nanmin(masked_new),
    "new_max": np.nanmax(masked_new),
    "new_mean": np.nanmean(masked_new),
    "new_median": np.nanmedian(masked_new),
    "new_std": np.nanstd(masked_new),
}
features

✔️ 131 ms (2024-12-02T11:04:15/2024-12-02T11:04:15)
Out[14]:
{'old_min': 264.148,
 'old_max': 520.245,
 'old_mean': 347.9752,
 'old_median': 350.69598,
 'old_std': 49.78903,
 'new_min': 264.706,
 'new_max': 519.89703,
 'new_mean': 356.3459,
 'new_median': 364.74097,
 'new_std': 53.15585}
In [15]:
raster = rasterize([largest_area.geometry], out_shape=result.shape, transform=new.transform)
masked_diff = np.where(raster, diff, np.nan)
masked_diff = masked_diff[top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_diff)

✔️ 1.82 s (2024-12-02T11:04:15/2024-12-02T11:04:17)
No description has been provided for this image
Out[15]:
<Axes: >
In [16]:
features.update({
    "diff_min": np.nanmin(masked_diff),
    "diff_max": np.nanmax(masked_diff),
    "diff_mean": np.nanmean(masked_diff),
    "diff_median": np.nanmedian(masked_diff),
    "diff_std": np.nanstd(masked_diff),
})
features

✔️ 78.1 ms (2024-12-02T11:04:17/2024-12-02T11:04:17)
Out[16]:
{'old_min': 264.148,
 'old_max': 520.245,
 'old_mean': 347.9752,
 'old_median': 350.69598,
 'old_std': 49.78903,
 'new_min': 264.706,
 'new_max': 519.89703,
 'new_mean': 356.3459,
 'new_median': 364.74097,
 'new_std': 53.15585,
 'diff_min': -6.8740234,
 'diff_max': 43.521973,
 'diff_mean': 8.370939,
 'diff_median': 4.911972,
 'diff_std': 9.06017}
In [17]:
attribute_names = ["roughness", "slope", "aspect", "curvature", "terrain_ruggedness_index", "rugosity", "profile_curvature", "planform_curvature"]
attributes = xdem.terrain.get_terrain_attribute(
    masked_diff,
    resolution=old.res,
    attribute=attribute_names
)

plt.figure(figsize=(8, 6.5))
for i in range(8):
    plt.subplot(4, 2, i + 1)
    plt.imshow(attributes[i].squeeze(), cmap="viridis")
    cbar = plt.colorbar()
    cbar.set_label(attribute_names[i])
    plt.xticks([])
    plt.yticks([])

plt.tight_layout()
plt.show()

✔️ 26.7 s (2024-12-02T11:04:17/2024-12-02T11:04:44)
No description has been provided for this image
In [18]:
for i, name in enumerate(attribute_names):
    features.update({
        f"{name}_min": np.nanmin(attributes[i]),
        f"{name}_max": np.nanmax(attributes[i]),
        f"{name}_mean": np.nanmean(attributes[i]),
        f"{name}_median": np.nanmedian(attributes[i]),
        f"{name}_std": np.nanstd(attributes[i]),
    })
features

✔️ 655 ms (2024-12-02T11:04:44/2024-12-02T11:04:45)
Out[18]:
{'old_min': 264.148,
 'old_max': 520.245,
 'old_mean': 347.9752,
 'old_median': 350.69598,
 'old_std': 49.78903,
 'new_min': 264.706,
 'new_max': 519.89703,
 'new_mean': 356.3459,
 'new_median': 364.74097,
 'new_std': 53.15585,
 'diff_min': -6.8740234,
 'diff_max': 43.521973,
 'diff_mean': 8.370939,
 'diff_median': 4.911972,
 'diff_std': 9.06017,
 'roughness_min': 0.024993896484375,
 'roughness_max': 16.99798583984375,
 'roughness_mean': 0.8664416971654487,
 'roughness_median': 0.501007080078125,
 'roughness_std': 0.9008248692669407,
 'slope_min': 0.00030909907704678905,
 'slope_max': 83.09449329265449,
 'slope_mean': 15.825287450127606,
 'slope_median': 9.967977228482425,
 'slope_std': 14.601294984565905,
 'aspect_min': 0.0,
 'aspect_max': 359.9997588896344,
 'aspect_mean': 172.74772778897523,
 'aspect_median': 183.00927436545985,
 'aspect_std': 104.09906838438863,
 'curvature_min': -1643.609619140625,
 'curvature_max': 2639.605712890625,
 'curvature_mean': 1.7315227064581844,
 'curvature_median': 0.8087158203125,
 'curvature_std': 43.72569863077058,
 'terrain_ruggedness_index_min': 0.024669455802839076,
 'terrain_ruggedness_index_max': 27.121265279417575,
 'terrain_ruggedness_index_mean': 0.8939862986187255,
 'terrain_ruggedness_index_median': 0.5406418612709043,
 'terrain_ruggedness_index_std': 0.8938326265050937,
 'rugosity_min': 1.000078699472661,
 'rugosity_max': 9.8867222334156,
 'rugosity_mean': 1.1083063933475281,
 'rugosity_median': 1.0261556432873709,
 'rugosity_std': 0.19958467669264315,
 'profile_curvature_min': -1767.2495203822261,
 'profile_curvature_max': 1433.1465250703304,
 'profile_curvature_mean': -1.2548662332856713,
 'profile_curvature_median': -0.5842965641336939,
 'profile_curvature_std': 28.49170269374419,
 'planform_curvature_min': -624.432616658446,
 'planform_curvature_max': 978.0291066330788,
 'planform_curvature_mean': 0.4766817029082818,
 'planform_curvature_median': 0.18651537440086147,
 'planform_curvature_std': 21.37279648106214}
In [19]:
REC = gpd.read_file("nzRec2_v5.gdb", engine='pyogrio', use_arrow=True)
REC

✔️ 3.4 s (2024-12-02T11:04:45/2024-12-02T11:04:48)
/home/ubuntu/.local/lib/python3.10/site-packages/pyogrio/raw.py:287: UserWarning: More than one layer found in 'nzRec2_v5.gdb': 'riverlines' (default), 'Hydro_Net_Junctions', 'rec2ws', 'rec2_rain_runoff_V5', 'rec2_stew_baserock_V5', 'rec2_stew_toprock_V5', 'rec2_utility_variables_V5', 'rec2_lcdb2_V5', 'rec2_lcdb1_V5', 'rec2_lcdb3_V5', 'rec2_ni_baserock_V5', 'rec2_ni_toprock_V5', 'rec2_si_baserock_V5', 'rec2_si_toprock_V5', 'rec2_elevband_rain_V5', 'LAYERKEYTABLE', 'APUNIQUEID', 'riverlines_FS', 'N_1_Desc', 'N_1_EDesc', 'N_1_EStatus', 'N_1_ETopo', 'N_1_FloDir', 'N_1_JDesc', 'N_1_JStatus', 'N_1_JTopo', 'N_1_JTopo2', 'N_1_Props'. Specify layer parameter to avoid this warning.
  with open_arrow(
Out[19]:
OBJECTID HydroID NextDownID CATAREA CUM_AREA nzsegment Enabled LENGTHDOWN Headwater Hydseq StreamOrde euclid_dis upElev downElev upcoordX downcoordX downcoordY upcoordY sinuosity nzreach_re headw_dist segslpmean LID reachtype FROM_NODE TO_NODE Shape_Leng Shape_Length FLOWDIR geometry
0 1 1 9 218553.786394 2.185538e+05 1000005 1 1801.260253 1 1 1 199.150697 119.994514 90.628487 1.601529e+06 1.601574e+06 6.192386e+06 6.192580e+06 1.072393 1000004 0 8.388180 0 0 1 2 213.567866 213.567866 1 MULTILINESTRING ((1601528.857 6192580.4, 16015...
1 2 2 9 455995.848576 4.559958e+05 1000003 1 1801.260253 1 2 1 537.313689 163.071075 90.628487 1.601783e+06 1.601574e+06 6.192386e+06 6.192881e+06 1.082775 1000005 0 7.678527 0 0 3 2 581.789877 581.789877 1 MULTILINESTRING ((1601782.856 6192881.076, 160...
2 3 3 4 461385.667778 4.613857e+05 1000008 1 335.656220 1 4 1 596.782205 64.632355 20.004650 1.599355e+06 1.599237e+06 6.191779e+06 6.192364e+06 1.163688 1000009 0 4.276652 0 0 4 5 694.468045 694.468045 1 MULTILINESTRING ((1599355.244 6192363.83, 1599...
3 4 4 -1 135806.633398 2.330307e+06 1000010 1 0.000000 0 9 2 288.064229 20.004650 4.172440 1.599237e+06 1.598982e+06 6.191913e+06 6.191779e+06 1.165213 1000010 788 3.145852 0 0 5 6 335.656220 335.656220 1 MULTILINESTRING ((1599237.074 6191778.666, 159...
4 5 5 -1 863421.445984 8.634214e+05 1000006 1 0.000000 1 5 1 1020.122542 143.468826 4.437851 1.602684e+06 1.602267e+06 6.191623e+06 6.192554e+06 1.182202 1000006 0 7.760943 0 0 7 8 1205.990464 1205.990464 1 MULTILINESTRING ((1602683.557 6192553.928, 160...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
593512 593614 593513 -1 51990.267885 1.174888e+09 4170001 1 NaN 0 124017 6 0.000000 0.000000 0.000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000 0 0 0.000000 0 0 85116 604841 542.278140 542.278140 1 MULTILINESTRING ((1975828.498 5785372.292, 197...
593513 593616 593514 -1 150666.238187 2.061735e+08 1040001 1 NaN 0 27152 5 0.000000 12.272364 0.000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000 0 0 0.000000 0 0 26864 604842 1675.015974 1675.015974 1 MULTILINESTRING ((1731848.082 6017486.385, 173...
593514 593620 593515 -1 53587.929573 5.503258e+07 7260002 1 NaN 0 239462 5 0.000000 4.917045 3.455193 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000 0 0 0.000000 0 0 244056 604843 156.868127 156.868127 1 MULTILINESTRING ((1789499.342 5528919.92, 1789...
593515 593621 593516 593515 8667.469549 5.432822e+07 7260001 1 NaN 0 239461 5 0.000000 6.187159 4.917045 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000 0 0 0.000000 0 0 244073 244056 116.148688 116.148688 1 MULTILINESTRING ((1789607.824 5528878.42, 1789...
593516 593622 593517 -1 492400.109562 1.856727e+08 7247115 1 529.317107 0 249508 5 964.116694 5.069054 1.886685 1.782739e+06 1.782650e+06 5.496834e+06 5.497794e+06 1.086063 7047051 36796 0.272352 0 0 252501 604844 151.702564 151.702564 1 MULTILINESTRING ((1782739.481 5497793.71, 1782...

593517 rows × 30 columns

In [20]:
LCDB = gpd.read_file("lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gpkg", engine='pyogrio', use_arrow=True)
LCDB

✔️ 6.03 s (2024-12-02T11:04:48/2024-12-02T11:04:54)
Out[20]:
Name_2018 Name_2012 Name_2008 Name_2001 Name_1996 Class_2018 Class_2012 Class_2008 Class_2001 Class_1996 Wetland_18 Wetland_12 Wetland_08 Wetland_01 Wetland_96 Onshore_18 Onshore_12 Onshore_08 Onshore_01 Onshore_96 EditAuthor EditDate LCDB_UID geometry
0 Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation 46 46 46 46 46 yes yes yes yes yes no no no no no Terralink 2004-06-30 lcdb2000096676 MULTIPOLYGON (((1613722.435 5425797.372, 16137...
1 Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation 46 46 46 46 46 yes yes yes yes yes no no no no no Regional Council 2019-12-01 lcdb1000513359 MULTIPOLYGON (((1816770.219 5947804.627, 18167...
2 Mangrove Mangrove Mangrove Mangrove Mangrove 70 70 70 70 70 yes yes yes yes yes no no no no no Terralink 2004-06-30 lcdb1000182160 MULTIPOLYGON (((1715672.186 5958842.706, 17156...
3 Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation Herbaceous Saline Vegetation 46 46 46 46 46 yes yes yes yes yes no no no no no Terralink 2004-06-30 lcdb1000065930 MULTIPOLYGON (((1705330.918 6088979.74, 170531...
4 Estuarine Open Water Estuarine Open Water Estuarine Open Water Estuarine Open Water Estuarine Open Water 22 22 22 22 22 no no no no no no no no no no Regional Council 2019-12-01 lcdb1000065472 MULTIPOLYGON (((1761684.636 5789742.527, 17616...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
511099 Low Producing Grassland Low Producing Grassland Low Producing Grassland Manuka and/or Kanuka Low Producing Grassland 41 41 41 52 41 no no no no no yes yes yes yes yes Landcare Research 2019-12-01 lcdb1000505810 MULTIPOLYGON (((1785112.194 5684560.595, 17851...
511100 Indigenous Forest Indigenous Forest Indigenous Forest Indigenous Forest Indigenous Forest 69 69 69 69 69 no no no no no yes yes yes yes yes Landcare Research 2019-12-01 lcdb2000193885 MULTIPOLYGON (((1607510.75 5432591.699, 160754...
511101 Exotic Forest Exotic Forest Exotic Forest Exotic Forest Exotic Forest 71 71 71 71 71 no no no no no yes yes yes yes yes Landcare Research 2011-06-30 lcdb2000219027 MULTIPOLYGON (((1603592.31 5269947.382, 160360...
511102 Herbaceous Freshwater Vegetation Herbaceous Freshwater Vegetation Herbaceous Freshwater Vegetation Herbaceous Freshwater Vegetation Herbaceous Freshwater Vegetation 45 45 45 45 45 yes yes yes yes yes yes yes yes yes yes Regional Council 2014-06-30 lcdb1000417326 MULTIPOLYGON (((1822629.596 5477414.336, 18226...
511103 High Producing Exotic Grassland High Producing Exotic Grassland Exotic Forest Exotic Forest Exotic Forest 40 40 71 71 71 no no no no no yes yes yes yes yes Landcare Research 2014-06-30 lcdb1000145507 MULTIPOLYGON (((1704402.741 5624819.729, 17044...

511104 rows × 24 columns

In [21]:
LCDB[LCDB.intersects(largest_area.geometry)]

✔️ 564 ms (2024-12-02T11:04:55/2024-12-02T11:04:55)
Out[21]:
Name_2018 Name_2012 Name_2008 Name_2001 Name_1996 Class_2018 Class_2012 Class_2008 Class_2001 Class_1996 Wetland_18 Wetland_12 Wetland_08 Wetland_01 Wetland_96 Onshore_18 Onshore_12 Onshore_08 Onshore_01 Onshore_96 EditAuthor EditDate LCDB_UID geometry
17220 Exotic Forest Exotic Forest Exotic Forest Exotic Forest High Producing Exotic Grassland 71 71 71 71 40 no no no no no yes yes yes yes yes Landcare Research 2011-06-30 lcdb1000053082 MULTIPOLYGON (((2054973.343 5805787.881, 20549...
18662 Gravel or Rock Gravel or Rock Gravel or Rock Gravel or Rock Gravel or Rock 16 16 16 16 16 no no no no no yes yes yes yes yes Landcare Research 2004-06-30 lcdb1000016122 MULTIPOLYGON (((2048258.705 5805949.487, 20482...
20295 Exotic Forest Exotic Forest Exotic Forest Exotic Forest Broadleaved Indigenous Hardwoods 71 71 71 71 54 no no no no no yes yes yes yes yes Landcare Research 2011-06-30 lcdb1000127897 MULTIPOLYGON (((2054496.624 5803468.778, 20545...
23069 Broadleaved Indigenous Hardwoods Broadleaved Indigenous Hardwoods Broadleaved Indigenous Hardwoods Broadleaved Indigenous Hardwoods Broadleaved Indigenous Hardwoods 54 54 54 54 54 no no no no no yes yes yes yes yes Landcare Research 2019-12-01 lcdb1000119141 MULTIPOLYGON (((2054140.516 5804625.89, 205416...
349500 Exotic Forest Exotic Forest Exotic Forest Exotic Forest Broadleaved Indigenous Hardwoods 71 71 71 71 54 no no no no no yes yes yes yes yes Landcare Research 2011-06-30 lcdb1000127894 MULTIPOLYGON (((2056301.039 5802553.702, 20562...
388204 Exotic Forest Exotic Forest Exotic Forest Exotic Forest Broadleaved Indigenous Hardwoods 71 71 71 71 54 no no no no no yes yes yes yes yes Landcare Research 2011-06-30 lcdb1000127900 MULTIPOLYGON (((2054727.599 5803954.788, 20547...
445455 Mixed Exotic Shrubland Mixed Exotic Shrubland Mixed Exotic Shrubland Mixed Exotic Shrubland Gravel or Rock 56 56 56 56 16 no no no no no yes yes yes yes yes Landcare Research 2019-12-01 lcdb2000515912 MULTIPOLYGON (((2055082.956 5803253.022, 20550...
In [30]:
LCDB[LCDB.intersects(largest_area.geometry)].drop(columns="EditDate").explore("Name_2018", legend=True, tiles="ESRI.WorldImagery")

✔️ 368 ms (2024-12-02T11:10:54/2024-12-02T11:10:54)
Out[30]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [23]:
LCDB.loc[LCDB.intersects(largest_area.geometry), ["Class_2018", "Wetland_18", "Onshore_18"]].mode().iloc[0].replace({"no": False, "yes": True})

✔️ 290 ms (2024-12-02T11:04:58/2024-12-02T11:04:58)
Out[23]:
Class_2018       71
Wetland_18    False
Onshore_18     True
Name: 0, dtype: object
In [24]:
all_streams = REC.unary_union

✔️ 10.4 s (2024-12-02T11:04:58/2024-12-02T11:05:08)
/tmp/ipykernel_3591039/3529544978.py:1: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.
  all_streams = REC.unary_union
In [25]:
features["distance_to_river"] = largest_area.geometry.distance(all_streams)
features["distance_to_river"]

✔️ 156 ms (2024-12-02T11:05:09/2024-12-02T11:05:09)
Out[25]:
0.0
In [26]:
def get_features(row):
    geom = row.geometry
    masked_old, transform = mask(old, [geom], nodata=np.nan)
    nonan = np.argwhere(~np.isnan(masked_old[0]))
    top_left = nonan.min(axis=0)
    bottom_right = nonan.max(axis=0)
    masked_old = masked_old[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]

    masked_new, transform = mask(new, [geom], nodata=np.nan)
    masked_new = masked_new[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]

    raster = rasterize([geom], out_shape=result.shape, transform=new.transform)
    masked_diff = np.where(raster, diff, np.nan)
    masked_diff = masked_diff[top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
    
    row = row.to_dict()
    row.update({
        "old_min": np.nanmin(masked_old),
        "old_max": np.nanmax(masked_old),
        "old_mean": np.nanmean(masked_old),
        "old_median": np.nanmedian(masked_old),
        "old_std": np.nanstd(masked_old),
        "new_min": np.nanmin(masked_new),
        "new_max": np.nanmax(masked_new),
        "new_mean": np.nanmean(masked_new),
        "new_median": np.nanmedian(masked_new),
        "new_std": np.nanstd(masked_new),
        "diff_min": np.nanmin(masked_diff),
        "diff_max": np.nanmax(masked_diff),
        "diff_mean": np.nanmean(masked_diff),
        "diff_median": np.nanmedian(masked_diff),
        "diff_std": np.nanstd(masked_diff),
        "distance_to_river": geom.distance(all_streams)
    })

    row.update(LCDB.loc[LCDB.intersects(geom), ["Class_2018", "Wetland_18", "Onshore_18"]].mode().iloc[0].replace({"no": False, "yes": True}))

    attribute_names = ["roughness", "slope", "aspect", "curvature", "terrain_ruggedness_index", "rugosity", "profile_curvature", "planform_curvature"]
    attributes = xdem.terrain.get_terrain_attribute(
        masked_diff,
        resolution=old.res,
        attribute=attribute_names
    )

    for i, name in enumerate(attribute_names):
        row.update({
            f"{name}_min": np.nanmin(attributes[i]),
            f"{name}_max": np.nanmax(attributes[i]),
            f"{name}_mean": np.nanmean(attributes[i]),
            f"{name}_median": np.nanmedian(attributes[i]),
            f"{name}_std": np.nanstd(attributes[i]),
        })

    return pd.Series(row)

#get_features(areas.iloc[0])
features = areas.progress_apply(get_features, axis=1)
features

✔️ 5 min 39 s (2024-12-02T11:05:09/2024-12-02T11:10:49)
  0%|          | 0/188 [00:00<?, ?it/s]
Out[26]:
geometry area old_min old_max old_mean old_median old_std new_min new_max new_mean new_median new_std diff_min diff_max diff_mean diff_median diff_std distance_to_river Class_2018 Wetland_18 Onshore_18 roughness_min roughness_max roughness_mean roughness_median roughness_std slope_min slope_max slope_mean slope_median slope_std aspect_min aspect_max aspect_mean aspect_median aspect_std curvature_min curvature_max curvature_mean curvature_median curvature_std terrain_ruggedness_index_min terrain_ruggedness_index_max terrain_ruggedness_index_mean terrain_ruggedness_index_median terrain_ruggedness_index_std rugosity_min rugosity_max rugosity_mean rugosity_median rugosity_std profile_curvature_min profile_curvature_max profile_curvature_mean profile_curvature_median profile_curvature_std planform_curvature_min planform_curvature_max planform_curvature_mean planform_curvature_median planform_curvature_std
176 POLYGON ((2054419 5803910, 2054419 5803909, 20... 530308.0 264.148010 520.244995 347.975189 350.695984 49.789028 264.705994 519.897034 356.345886 364.740967 53.155849 -6.874023 43.521973 8.370939 4.911972 9.060170 0.000000 71 False True 0.024994 16.997986 0.866442 0.501007 0.900825 0.000309 83.094493 15.825287 9.967977 14.601295 0.000000 359.999759 172.747728 183.009274 104.099068 -1643.609619 2639.605713 1.731523 0.808716 43.725699 0.024669 27.121265 0.893986 0.540642 0.893833 1.000079 9.886722 1.108306 1.026156 0.199585 -1767.249520 1433.146525 -1.254866 -0.584297 28.491703 -624.432617 978.029107 0.476682 0.186515 21.372796
155 POLYGON ((2054514 5804756, 2054514 5804755, 20... 319463.0 384.855011 949.697021 698.896851 715.611023 132.071274 384.348999 949.190002 684.987671 692.672974 129.693771 -52.056030 5.609009 -13.909487 -9.703979 12.755947 0.000000 54 False True 0.027039 10.799988 1.465731 1.246033 0.944709 0.021028 77.297658 26.350639 24.870094 14.383982 0.000000 359.999410 183.019811 185.773786 100.256831 -954.620361 692.401123 -1.137023 -1.690674 58.625323 0.035931 10.752627 1.490565 1.282474 0.923966 1.000176 4.733187 1.209573 1.132283 0.229804 -440.434936 585.141157 0.664059 0.707163 35.661075 -369.479204 312.140260 -0.473099 -0.604222 31.241830
37 POLYGON ((2056589 5808659, 2056589 5808657, 20... 304243.0 543.239014 968.459961 752.026978 758.137024 89.943939 542.731018 968.713989 750.921082 757.086975 90.210747 -13.114990 5.141968 -1.105809 -0.944031 1.355662 0.000000 12 False True 0.014038 6.578003 0.709351 0.576965 0.492890 0.019821 66.326472 13.694865 11.374013 9.495378 0.000000 359.998682 180.263642 180.474840 104.744639 -472.998047 416.595459 -1.971487 -1.379395 39.778407 0.014390 5.968466 0.742228 0.609774 0.502908 1.000017 2.568007 1.063572 1.033145 0.087242 -275.317190 316.405986 1.262352 0.660475 23.922942 -335.581215 226.580603 -0.709152 -0.401965 21.456914
115 POLYGON ((2054495 5806342, 2054495 5806340, 20... 173781.0 333.551025 461.656006 378.365387 375.186005 26.501810 334.056000 461.194000 381.083679 378.898010 26.786768 -4.665009 12.177979 2.718257 2.453003 1.906919 0.000000 12 False True 0.028015 6.747986 0.664326 0.454987 0.605888 0.028203 69.110675 12.514475 8.870377 11.030167 0.000000 359.999000 164.568461 177.295377 100.462405 -553.802490 687.600708 2.607275 1.797485 41.001715 0.032552 7.412247 0.699980 0.489235 0.623502 1.000124 2.874868 1.065598 1.021682 0.119312 -479.441262 366.685717 -2.050307 -1.421865 27.232389 -357.169244 403.752029 0.556968 0.261994 19.498026
126 POLYGON ((2054821 5806083, 2054821 5806082, 20... 147413.0 405.266998 984.378967 666.230042 680.250000 128.762848 404.644012 984.406982 662.956543 677.700012 128.024765 -25.545959 10.151001 -3.273446 -2.276978 3.636282 0.000000 71 False True 0.049988 8.898010 1.275577 1.083984 0.813165 0.040820 73.866107 23.143639 21.255515 13.418823 0.000000 359.998574 184.088882 176.004473 102.389142 -767.303467 707.580566 -4.001010 -3.497314 62.809494 0.051240 8.480720 1.322321 1.132742 0.821058 1.000402 3.628650 1.171587 1.106380 0.191066 -536.233840 418.001613 2.441112 1.776779 38.928714 -393.931141 379.570854 -1.559898 -1.174689 33.025140
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
51 POLYGON ((2057180 5807155, 2057180 5807154, 20... 4309.0 754.164001 862.838013 810.229065 812.798035 27.570267 753.013000 862.265991 808.988464 811.622009 27.611164 -5.234009 0.836975 -1.240615 -1.079041 0.706411 377.665475 69 False True 0.091980 3.499023 0.810893 0.724518 0.466198 0.258646 54.086741 15.535889 14.296695 9.054989 0.000000 359.873676 180.929910 181.396905 105.351676 -295.800781 442.993164 -7.090177 -4.492188 47.193806 0.086022 3.618491 0.855655 0.758100 0.486536 1.001061 1.724303 1.077381 1.049964 0.088696 -238.362703 175.115296 2.803598 1.571423 27.481605 -176.705037 204.630461 -4.286579 -2.724701 26.676936
77 POLYGON ((2057221 5806249, 2057221 5806248, 20... 4231.0 463.237000 558.276001 506.580353 501.701996 21.162735 462.401001 557.640991 504.420471 500.536987 21.546041 -9.316010 0.420044 -2.159984 -1.524994 1.754438 0.000000 16 False True 0.123962 5.431030 1.354514 1.208984 0.781550 0.243961 63.330647 24.899263 23.870917 12.951836 0.117864 359.932983 184.918589 175.719571 102.218871 -365.005493 229.119873 -11.565558 -9.808350 54.859841 0.102363 4.997128 1.382453 1.254762 0.757155 1.000989 2.266732 1.179069 1.126430 0.172435 -180.493919 211.376430 5.078110 3.214066 31.870435 -207.633207 137.569661 -6.487448 -5.183216 31.373435
181 POLYGON ((2054921 5802132, 2054921 5802131, 20... 4168.0 727.848999 777.617981 750.956482 751.011475 10.164408 727.346008 777.073975 748.487061 747.432495 10.804439 -7.746033 0.098999 -2.469438 -1.833984 1.776451 313.839735 71 False True 0.050049 2.488037 0.713715 0.641541 0.372780 0.161113 41.709637 14.652946 13.619203 7.531259 0.239895 359.927197 183.505383 172.413544 104.640565 -209.301758 189.300537 -3.036568 -1.809692 25.295615 0.050731 2.201049 0.710312 0.651643 0.352425 1.000237 1.344907 1.051155 1.036017 0.048002 -116.756886 137.333515 1.427443 0.626500 14.888978 -90.982558 93.936800 -1.609124 -0.725124 13.837087
112 POLYGON ((2056165 5805098, 2056165 5805097, 20... 4168.0 330.009003 414.976013 347.387390 338.135986 19.942204 329.310974 414.473022 345.676178 335.724487 20.337126 -6.242035 0.404999 -1.711227 -1.103989 1.331259 20.061887 16 False True 0.067993 5.648010 0.916036 0.627975 0.812676 0.072869 67.422201 16.667242 11.679819 14.067467 0.000371 359.989076 160.836444 161.757142 101.798624 -645.401001 310.797119 -10.793974 -5.749512 54.868025 0.061296 7.057928 0.977414 0.664636 0.868992 1.000381 2.796020 1.118655 1.040430 0.184431 -213.102706 411.064214 7.399746 3.237263 38.168409 -249.806248 200.170430 -3.394228 -2.160090 26.729371
70 POLYGON ((2056954 5806554, 2056954 5806551, 20... 4069.0 580.815002 632.517029 607.255920 607.020020 11.594285 580.307007 632.016968 605.696106 604.859009 11.600400 -5.670044 -0.108032 -1.559806 -1.331970 0.911266 44.852111 69 False True 0.078003 2.404053 0.727969 0.664001 0.367533 0.283751 40.625595 14.404642 13.091004 7.748489 0.185128 359.981782 184.885313 190.514454 104.413637 -207.397461 179.193115 -6.322378 -4.388428 34.206604 0.078687 2.184222 0.748892 0.687885 0.359458 1.000577 1.326722 1.056693 1.041069 0.049973 -87.630775 112.695939 2.519276 1.309808 19.584720 -115.848622 91.562340 -3.803102 -2.074324 19.314333

188 rows × 61 columns

In [27]:
features.crs = new.crs

✔️ 157 ms (2024-12-02T11:10:49/2024-12-02T11:10:50)
In [28]:
features.explore("Class_2018", legend=True, tiles="ESRI.WorldImagery")

✔️ 1.95 s (2024-12-02T11:10:50/2024-12-02T11:10:52)
Out[28]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [29]:
features.to_parquet("BD44_10000_0503.parquet")

✔️ 64.3 ms (2024-12-02T11:10:54/2024-12-02T11:10:54)